indTaRKD
Quality assessment of trimmed reads
FastQC
The reads were trimmed using fastp an ultra-fast all-in-one FASTQ preprocessor" and quality checked using FASTQC. All quality report plots were collated using MultiQC and plots were exported and visualised below. No irregularities were found in the trimmed reads so we continued to alignment.
Figure 1A. FastQC per sequence quality score plot
Figure 1B. FastQC per sequence quality score plot
Figure 1C. FastQC per sequence quality score plot
Figure 1D. FastQC per sequence quality score plot
Kallisto pseudoalignment
Trimmed reads were mapped to the RefSeq v1.1 wheat genome using a pseudoaligner called kallisto. Kallisto is much faster than conventional mapping tools and has been shown to have equal or greater mapping efficiency because it skips the alignment step (exact position) and just focuses on mapping to the best match in a indexed cDNA library. Below is a graphic of the number of reads processed/aligned successfully.
Figure 1E. Read_processed plot
Figure 1F. Percentage pseudoaligned plot
DESeq2 normalisation and analysis of replicate quality
indTaRKD dispersion plot - DESeq2 normalisation
Counts were loaded into R and normalised using DESeq2. Below is a dispersion plot where we can see that counts were fitted to the trend of gene expression dispersion (variance). Outliers (bottom of the plot) are excluded from the analysis.
Figure 2A. Wheat dispersion plot
Principal component analysis of normalised counts
Clear seperation of transgenic and non-transgenic samples along PC1 which accounts for 83% of gene expression dispersion (variance).
Figure 2B. indTaRKD PCA
Scatter plots of replicates within conditions
Scatterplots of normalised counts for genes. We can see good concordance between replicates suggesting good biological replication.
Figure 2C. indTaRKD replicate scatterplots
DESeq2 analysis of expression profiles
indTaRKD expression
We could see a clear difference in indTaRKD expression across samples that we hypothesise is causative of the tillering phenotype.
Figure 3A. Normalised indTaRKD expression across sample
Diffrential expression analysis
DESeq2 calculated diffrentially expressed genes under the following criteria; lfcThreshold = 0.5, p.value = 0.05.
## out.of.2979.with.nonzero.total.read.count
## 1 adjusted p-value < 0.05
## 2 LFC > 0.50 (up) : 1494, 50%
## 3 LFC < -0.50 (down) : 1485, 50%
## 4 outliers [1] : 0, 0%
## 5 low counts [2] : 0, 0%
## 6 (mean count < 5)MA plots to see pattern of up/down regulation of normalised counts
Equal amount of up and downregulation, consistent with this being
Figure 3C. indTaRKD MA plot Transgenic vs Non-transgenic
Heatmap of diffrentially expressed library reads
Heatmap of logarithmically transformed normalised counts of indTaRKD libraries. We can see clear differences in upregulated and downregulated genes between samples
Figure 3D. indTaRKD DEGs Heatmaps
Enrichment of functional terms
Transcription factor enrichment
Using a package called gage we tested for transcription factor enrichment. We used an enrichment strategy that takes normalised counts as an input as opposed to DEGs. This allows us to see more accurate enrichment of classes of transcription factors. One limitation of this strategy is that large transcription factor familes will not be enriched (plot below).
Figure 4A. Set size vs. enrichment rank
## p.val q.val set.size
## C3H 1.273463e-15 7.258737e-14 89
## TUB 6.650340e-13 1.895347e-11 32
## Sigma70-like 2.509744e-10 4.768514e-09 13
## TAZ 1.451101e-09 2.067819e-08 16
## CAMTA 2.078538e-09 2.369534e-08 13
## Pseudo ARR-B 2.846642e-09 2.704310e-08 14
## ARF 6.898950e-09 5.617716e-08 50
## C2C2_CO-like 1.574470e-08 1.121810e-07 22
## tify 3.055970e-06 1.935448e-05 54
## HSF 4.360328e-06 2.485387e-05 60
## HD_PHD 5.037393e-06 2.610285e-05 9
## ARID 2.680387e-05 1.217615e-04 9
## Zinc finger, MIZ type 2.777017e-05 1.217615e-04 12
## E2F/DP 7.560278e-05 3.078113e-04 23
## EIL 1.117401e-04 4.246123e-04 14
## FHA 1.411257e-04 5.027602e-04 41
## VOZ 2.205378e-04 7.394502e-04 5
## GeBP 2.984212e-04 9.450004e-04 33
## RB 4.570084e-04 1.371025e-03 5
## BES1 1.890849e-03 5.388919e-03 16
## Sir2 1.511833e-02 4.103547e-02 4
## Whirly 1.659774e-02 4.300324e-02 5
## HRT 3.203280e-02 7.938564e-02 2
## RWP-RK 3.550248e-02 8.431838e-02 27
## CCAAT_HAP2 4.247663e-02 9.684671e-02 14
## CSD 5.273003e-02 1.156005e-01 12
## BBR/BPC 1.700229e-01 3.589371e-01 4
## ULT 2.730318e-01 5.372288e-01 3
## C2C2_GATA 2.733269e-01 5.372288e-01 64
## WRKY 3.271276e-01 6.215425e-01 244
## CPP 3.876729e-01 6.950320e-01 32
## bHLH_TCP 3.901934e-01 6.950320e-01 44
## GARP_G2-like 4.628251e-01 7.994252e-01 84
## GRAS 5.058029e-01 8.479637e-01 145
## CCAAT_HAP3 5.243669e-01 8.539690e-01 2
## HD-Zip_I_II 5.736520e-01 9.082823e-01 43
## CCAAT_HAP5 6.425949e-01 9.899435e-01 29
## HD-Zip_IV 6.614274e-01 9.921411e-01 33
## MYB-related 6.950654e-01 1.000000e+00 506
## C2C2_Dof 7.101377e-01 1.000000e+00 82
## HD 7.995020e-01 1.000000e+00 170
## CCAAT_Dr1 8.058911e-01 1.000000e+00 30
## mTERF 8.872037e-01 1.000000e+00 229
## HD_PLINC 8.884779e-01 1.000000e+00 32
## SBP 9.652327e-01 1.000000e+00 44
## SRS 9.800760e-01 1.000000e+00 13
## LFY 9.891216e-01 1.000000e+00 3
## C2C2_YABBY 9.921706e-01 1.000000e+00 16
## GRF 9.978584e-01 1.000000e+00 27
## PLATZ 9.994254e-01 1.000000e+00 39
## AS2/LOB 9.999992e-01 1.000000e+00 68
## AP2/EREBP 9.999998e-01 1.000000e+00 405
## bHLH 1.000000e+00 1.000000e+00 277
## MADS_II 1.000000e+00 1.000000e+00 108
## NAC 1.000000e+00 1.000000e+00 347
## ABI3/VP1 1.000000e+00 1.000000e+00 244
## MADS_I 1.000000e+00 1.000000e+00 233
## S1Fa-like NA NA 1Signficantly enriched transcription factors include processes like signalling and polymerase activity like CAMTA and Sigma70. Sigma70 TFs are involved in RNA polymerisation and are highly significanly upregulated in our indTaRKD DEGs. Additionally they’re upregulated in our OsRKD and microspore datasets and have a predicted RKD motif in their promoter region. Other TFs that are consistent with this is E2F-DP which are involved in the G1/S transition.
C3H transcription factors are involved in photomorphogenesis early in embryogenesis. Consistent with this photosynthesis role is CO-like TFs have also been shown to be involved in flowering time.
The final point is that BES1 and tify transcription factors are involved in brassinosteroid and jasmonate response pathways that have been shown to be involved in tillering and regeneration.
All of this is suggests a transcriptional cascade that regulateds cell division and transcription which requires energy from photosynthesis. Two phytohormones involved in the tillering and regeneration context of RKD is Brassinosteroids and Jasmonate go onto explain our phenotype.
Figure 4C. Normalised TaDWF4 (BR synthesis) expression across samples.
Go term enrichment
Over-representation
To conduct a GO term analysis we took the RefSeq v1.0 GO terms list made available by Philippa Borrill and translated these to RefSeq v1.1 IDs (99% kept). Then we used the package TopGO to take compare GO term observed in our DEGs with the expected number.
The first plot shows enrichment of GO terms involved in energy production and cell division. Xanathine is invovled in purine metabolism and shows very high enrichment along with terms involved in sucrose transport and starch catabolism. Other interesting terms include brassinosteroid biosynthesis, response to wounding, defense response, cytokinin metabolic proccesses along with many photosynthesis related terms.
Figure 4D. indTaRKD Top 50 enriched GO term wheat ratio
Under-representation
The second plot shows depletion of GO-terms involved in ubiquitin mediated degredation, regulation of transcription, innate immune response and phosphorylation. Can’t see any obvious trends from this other than perhaps a diversion of resources from immune related process to growth and development. Under-representation is usually not very relevant however…
Figure 4E. indTaRKD Top 50 depleted GO term wheat ratio
The two plots below GO term analysis using -log10(p-values). This gives a better idea of the most significant terms that are enriched/depleted
Figure 4F. indTaRKD Top40 enriched GO terms wheat -log10(p-value)
Figure 4G. indTaRKD Top 40 depleted GO terms wheat -log10(p-value)
Motif enrichment of 500bp upstream promoter sequences of upregulated DEGs
500bp promoter sequences were extracted from the RefSeq v1.1 assembly using the R package GRanges. Upregulated DEG sequences were compared against a background group of randomly chosen genes and analysed using HOMER motif enrichment analysis using the findMotifs.pl fasta command. This generated a html output of the de-novo and known motifs detected.
Figure 4H. HOMER Known motif enrichment in 500bp upstream indTaRKD upregulated genes
We see high enrichment of a MYB transcription factor motif. These TFs are usually involved in regulating cell division and diffrentiation. A high confidence candidate in our indTaRKD DEGs is a MYB TF that is upregulated in a microspore dataset we’ve analysed (described later) and has a predicted RKD motif in its promoter region. The 5th motif is a Phytoclock1 related motif consistent with many downstream targets having roles in photomorphogenesis.
Figure 4I. Normalised MYB10 expression across samples
Figure 4J. Normalised PCL1 triad expression across samples
Figure 4K. HOMER de-novo motif enrichment in 500bp upstream indTaRKD upregulated genes
The only de-novo motif that crosses the e-12 threshold (the cut-off for not risking false positives) is an AP2-ERF motif. These transcription factors have been implicated in regeneration and response to abiotic stress. There is an AP2-ERF significantly upregulated in the both the OsRKD and indTaRKD DEGs (ERF113) that has an RKD-motif in the rice homolog. Additionally this AP2-ERF has been shown to play a role in wound response and regeneration in the literature.
Figure 4L. Normalised expression of ERF113 (DEG) and its homeolog ERF115 across samples (non-DEG)
Microspore
Experimental design
This study took the bulgarian winter wheat ‘Svilena’ which is known for its high regenerative properties and took three tissue samples; untreated microspores (non-dividing), pretreated (cold stress - non-dividing) and treated (dividing) microspores. We decided to analyse this dataset as it had not been done with the RefSeq assembly and it would provide a novel insight into early embryogenesis in wheat.
DESeq2 normalisation and analysis of replicate quality
indTaRKD dispersion plot - DESeq2 normalisation
Counts were loaded into R and normalised using DESeq2. Below is a dispersion plot where we can see that counts were fitted to the trend of gene expression dispersion (variance). Outliers (bottom of the plot) are excluded from the analysis.
Figure 5A. Microspore dispersion plot
Principal component analysis of normalised counts
Clear seperation of transgenic and non-transgenic samples along PC1 which accounts for 94% of gene expression dispersion (variance).
Figure 5B. Microspore PCA
Diffrential expression analysis
DESeq2 calculated diffrentially expressed genes under the following criteria; lfcThreshold = 2, p.value = 0.01.
## out.of.85066.with.nonzero.total.read.count
## 1 adjusted p-value < 0.01
## 2 LFC > 2.00 (up) : 4349, 5.1%
## 3 LFC < -2.00 (down) : 2384, 2.8%
## 4 outliers [1] : 395, 0.46%
## 5 low counts [2] : 32465, 38%
## 6 (mean count < 8)## out.of.85066.with.nonzero.total.read.count
## 1 adjusted p-value < 0.01
## 2 LFC > 2.00 (up) : 3313, 3.9%
## 3 LFC < -2.00 (down) : 1768, 2.1%
## 4 outliers [1] : 395, 0.46%
## 5 low counts [2] : 34047, 40%
## 6 (mean count < 9)## out.of.5081.with.nonzero.total.read.count
## 1 adjusted p-value < 0.01
## 2 LFC > 2.00 (up) : 2937, 58%
## 3 LFC < -2.00 (down) : 1353, 27%
## 4 outliers [1] : 0, 0%
## 5 low counts [2] : 0, 0%
## 6 (mean count < 8)Venn diagram of OsRKD DEGs
Equal amount of up and downregulation, consistent with this being
Figure 5F. Venn diagram of Microspore DEGs
Heatmap of diffrentially expressed library reads
Heatmap of logarithmically transformed normalised counts of microspore libraries. We can see clear differences in upregulated and downregulated genes between indOsRKD expressing and non-expressing samples
Figure 5G. Shared Microspore DEGs Heatmaps
Enrichment of functional terms
Transcription factor enrichment
Using a package called gage we tested for transcription factor enrichment. We used an enrichment strategy that takes normalised counts as an input as opposed to DEGs. This allows us to see more accurate enrichment of classes of transcription factors. One limitation of this strategy is that large transcription factor familes will not be enriched (plot below).
## p.val q.val set.size
## ARF 2.457723e-16 1.400902e-14 68
## C3H 8.759784e-15 2.496538e-13 101
## CAMTA 3.934844e-11 7.476204e-10 16
## ARID 1.101425e-07 1.569531e-06 13
## Pseudo ARR-B 4.635353e-07 5.284303e-06 14
## Zinc finger, MIZ type 7.337862e-07 6.970969e-06 15
## FHA 7.817136e-06 6.365382e-05 49
## HD_PHD 1.167263e-05 7.880100e-05 9
## E2F/DP 1.244226e-05 7.880100e-05 27
## TUB 1.670110e-04 9.519626e-04 36
## Whirly 8.702589e-04 4.509523e-03 6
## CPP 1.678203e-03 7.955558e-03 35
## VOZ 1.814425e-03 7.955558e-03 6
## S1Fa-like 5.698827e-03 2.320237e-02 3
## RB 3.253437e-02 1.236306e-01 6
## Sir2 4.131658e-02 1.471903e-01 5
## TAZ 8.987753e-02 3.013541e-01 20
## HRT 1.062143e-01 3.363452e-01 2
## CCAAT_HAP3 1.282799e-01 3.848398e-01 3
## EIL 1.759560e-01 4.970889e-01 14
## GeBP 1.831380e-01 4.970889e-01 36
## Sigma70-like 2.210594e-01 5.727447e-01 18
## CCAAT_HAP2 2.656035e-01 6.582347e-01 15
## BBR/BPC 3.001473e-01 7.128498e-01 6
## CCAAT_Dr1 3.315523e-01 7.559392e-01 17
## MADS_II 4.445534e-01 9.745978e-01 108
## SBP 4.824579e-01 1.000000e+00 45
## RWP-RK 5.161537e-01 1.000000e+00 37
## tify 7.175362e-01 1.000000e+00 56
## HSF 7.356910e-01 1.000000e+00 58
## BES1 8.262842e-01 1.000000e+00 15
## GRF 8.925792e-01 1.000000e+00 21
## CCAAT_HAP5 9.362726e-01 1.000000e+00 28
## ULT 9.670014e-01 1.000000e+00 3
## C2C2_YABBY 9.716900e-01 1.000000e+00 21
## C2C2_CO-like 9.965979e-01 1.000000e+00 13
## CSD 9.970570e-01 1.000000e+00 15
## mTERF 9.992952e-01 1.000000e+00 294
## HD-Zip_IV 9.994769e-01 1.000000e+00 29
## C2C2_GATA 9.999859e-01 1.000000e+00 55
## SRS 9.999965e-01 1.000000e+00 9
## ABI3/VP1 9.999985e-01 1.000000e+00 251
## C2C2_Dof 9.999990e-01 1.000000e+00 48
## GRAS 9.999998e-01 1.000000e+00 128
## MADS_I 1.000000e+00 1.000000e+00 185
## MYB-related 1.000000e+00 1.000000e+00 449
## AS2/LOB 1.000000e+00 1.000000e+00 43
## GARP_G2-like 1.000000e+00 1.000000e+00 76
## HD 1.000000e+00 1.000000e+00 147
## HD_PLINC 1.000000e+00 1.000000e+00 22
## HD-Zip_I_II 1.000000e+00 1.000000e+00 39
## NAC 1.000000e+00 1.000000e+00 317
## bHLH_TCP 1.000000e+00 1.000000e+00 45
## PLATZ 1.000000e+00 1.000000e+00 24
## AP2/EREBP 1.000000e+00 1.000000e+00 277
## bHLH 1.000000e+00 1.000000e+00 213
## WRKY 1.000000e+00 1.000000e+00 171
## LFY NA NA 0
Figure 6A. indTaRKD Transcription factor enrichment
Similar TF enrichment to indTaRKD with C3H, CAMTA, E2F, and TUB playing roles. The big noteable difference being ARF transcription factors which is consistent with microspore induction (the first division) being highly dependent on auxin.
Go term enrichment
Over-representation
To conduct a GO term analysis we took the RefSeq v1.0 GO terms list made available by Philippa Borrill and translated these to RefSeq v1.1 IDs (99% kept). Then we used the package TopGO to take compare GO term observed in our DEGs with the expected number.
Jasmonate response appears in enriched GO terms but brassinosteroid does not. This suggests that brassinosteroid isn’t the target of RKD (a zygotic factor) but that jasmonate could be.
Figure 6B. Microspore Top 50 enriched GO term wheat ratio
Figure 6C. Microspore Top 50 depleted GO term wheat ratio
Motif enrichment of 500bp upstream promoter sequences of upregulated DEGs
500bp promoter sequences were extracted from the RefSeq v1.1 assembly using the R package GRanges. Upregulated DEG sequences were compared against a background group of randomly chosen genes and analysed using HOMER motif enrichment analysis using the findMotifs.pl fasta command. This generated a html output of the de-novo and known motifs detected.
Figure 6D. HOMER Known motif enrichment in 500bp upstream indTaRKD upregulated genes
There is a high enrichment of MYB transcription factor motifs and no de-novo enriched motifs.
indOsRKD
DESeq2 normalisation and analysis of replicate quality
indOsRKD dispersion plot - DESeq2 normalisation
Counts were loaded into R and normalised using DESeq2. Below is a dispersion plot where we can see that counts were fitted to the trend of gene expression dispersion (variance). Outliers (bottom of the plot) are excluded from the analysis.
Figure 7A. indOsRKD dispersion plot
Principal component analysis of normalised counts
Clear seperation of OsRKD induced samples along PC1 which accounts for 92% of gene expression dispersion (variance). Slight seperation of samples according to genotype along PC2 (~2%).
Figure 7B. indTaRKD PCA
DESeq2 analysis of expression profiles
indOsRKD expression
We could see a clear difference in indTaRKD expression across samples that we hypothesise is causative of the tillering phenotype.
Figure 8A. Normalised indTaRKD expression across sample
Diffrential expression analysis
DESeq2 calculated diffrentially expressed genes under the following criteria; lfcThreshold = 1, p.value = 0.05.
## out.of.4280.with.nonzero.total.read.count
## 1 adjusted p-value < 0.05
## 2 LFC > 1.00 (up) : 1640, 38%
## 3 LFC < -1.00 (down) : 2640, 62%
## 4 outliers [1] : 0, 0%
## 5 low counts [2] : 0, 0%
## 6 (mean count < 4)## out.of.3928.with.nonzero.total.read.count
## 1 adjusted p-value < 0.05
## 2 LFC > 1.00 (up) : 1622, 41%
## 3 LFC < -1.00 (down) : 2306, 59%
## 4 outliers [1] : 0, 0%
## 5 low counts [2] : 0, 0%
## 6 (mean count < 4)## out.of.3259.with.nonzero.total.read.count
## 1 adjusted p-value < 0.05
## 2 LFC > 1.00 (up) : 1210, 37%
## 3 LFC < -1.00 (down) : 2049, 63%
## 4 outliers [1] : 0, 0%
## 5 low counts [2] : 0, 0%
## 6 (mean count < 4)Venn diagram of OsRKD DEGs
Equal amount of up and downregulation, consistent with this being
Figure 8E. Venn diagram of OsRKD DEGs
Heatmap of diffrentially expressed library reads
Heatmap of logarithmically transformed normalised counts of indOsRKD libraries. We can see clear differences in upregulated and downregulated genes between indOsRKD expressing and non-expressing samples
Figure 8F. indOsRKD DEGs Heatmaps
Enrichment of functional terms
Transcription factor enrichment
Using a package called gage we tested for transcription factor enrichment. We used an enrichment strategy that takes normalised counts as an input as opposed to DEGs. This allows us to see more accurate enrichment of classes of transcription factors. One limitation of this strategy is that large transcription factor familes will not be enriched (plot below).
## p.val q.val set.size
## CO-like 2.274682e-07 1.228329e-05 13
## WRKY 1.798076e-05 4.854806e-04 51
## CAMTA 4.438055e-04 7.988499e-03 5
## DBB 1.001489e-03 1.352010e-02 7
## HB-other 5.290208e-03 5.713425e-02 5
## RAV 1.603369e-02 1.443032e-01 3
## MYB_related 2.514382e-02 1.939666e-01 24
## VOZ 3.819955e-02 2.578470e-01 2
## NF-X1 4.841947e-02 2.905168e-01 2
## ARR-B 5.646091e-02 3.048889e-01 4
## ARF 6.868628e-02 3.371872e-01 12
## G2-like 9.123926e-02 4.105767e-01 29
## Nin-like 1.083592e-01 4.501072e-01 7
## CPP 1.292174e-01 4.792709e-01 3
## BES1 1.331308e-01 4.792709e-01 5
## Trihelix 1.574149e-01 5.312753e-01 19
## C3H 1.838222e-01 5.839057e-01 24
## S1Fa-like 2.149447e-01 6.308513e-01 2
## Whirly 2.407433e-01 6.308513e-01 2
## EIL 2.408574e-01 6.308513e-01 3
## LSD 2.506702e-01 6.308513e-01 5
## GRAS 2.570135e-01 6.308513e-01 28
## Dof 4.147610e-01 9.737867e-01 18
## NF-YB 4.594658e-01 1.000000e+00 7
## ERF 5.471179e-01 1.000000e+00 78
## GATA 6.910952e-01 1.000000e+00 14
## HSF 7.262724e-01 1.000000e+00 15
## NF-YA 7.413524e-01 1.000000e+00 2
## TALE 7.857473e-01 1.000000e+00 13
## NF-YC 8.126593e-01 1.000000e+00 8
## E2F/DP 8.646824e-01 1.000000e+00 6
## SBP 8.650753e-01 1.000000e+00 12
## NAC 8.770866e-01 1.000000e+00 62
## FAR1 8.932826e-01 1.000000e+00 23
## ZF-HD 9.162308e-01 1.000000e+00 14
## TCP 9.304631e-01 1.000000e+00 13
## HD-ZIP 9.334414e-01 1.000000e+00 24
## BBR-BPC 9.407126e-01 1.000000e+00 2
## GeBP 9.681360e-01 1.000000e+00 13
## bZIP 9.694513e-01 1.000000e+00 46
## GRF 9.711163e-01 1.000000e+00 5
## YABBY 9.831112e-01 1.000000e+00 3
## M-type_MADS 9.858032e-01 1.000000e+00 5
## bHLH 9.969956e-01 1.000000e+00 73
## WOX 9.982648e-01 1.000000e+00 3
## MIKC_MADS 9.995228e-01 1.000000e+00 16
## SRS 9.997595e-01 1.000000e+00 4
## MYB 9.999945e-01 1.000000e+00 111
## AP2 9.999990e-01 1.000000e+00 7
## C2H2 1.000000e+00 1.000000e+00 63
## LBD 1.000000e+00 1.000000e+00 21
## B3 1.000000e+00 1.000000e+00 25
## LFY NA NA 1
## HB-PHD NA NA 1
## HRT-like NA NA 1Go term enrichment
To conduct a GO term analysis we took GoSlim terms from the MSU database and plotted significant values similar to before.
Over-representation
Figure 9B. Top 10 enriched GOSlim rice ratio
We see a high enrichment of GOslim terms related to photosynthesis, energy related metabolic processes and response to abiotic stress.
Under-representation
Figure 9C. Top 10 depleted GOSlim rice ratio
Motif enrichment of 1000bp upstream promoter sequences of upregulated DEGs
1000bp upstream promoter sequences were taken from the RAP-DB database. Upregulated Shared DEG (p.value < 0.01, lfcThreshold > 1 - 618 DEGs) sequences were compared against a background group of randomly chosen genes and analysed using HOMER motif enrichment analysis using the findMotifs.pl fasta command. This generated a html output of the de-novo and known motifs detected.
Figure 9D. HOMER Known motif enrichment in 500bp upstream indOsRKD upregulated genes
The only de-novo motif that crosses the e-12 threshold (the cut-off for not risking false positives) is an AP2-ERF motif. These transcription factors have been implicated in regeneration and response to abiotic stress. There is an AP2-ERF significantly upregulated in the both the OsRKD and indTaRKD DEGs (ERF113) that has an RKD-motif in the rice homolog. Its expression is brassinosteroid induced. Additionally this AP2-ERF has been shown to play a role in wound response and regeneration in the literature.
Figure 9E. HOMER Known motif enrichment in 500bp upstream indOsRKD upregulated genes
Figure 9F. HOMER Known motif enrichment in 500bp upstream indOsRKD upregulated genes
In the known motifs there is a high enrichment of WRKY, MYB and AP2-ERF transcription factors. Of particular interest is the enrichment of ERF115, a homolog to the ERF113 homeolog that controls root QC cell division and stem cell niche replenishment. All of this once again points to a role of ERFs regulating regeneration and meristem initiation/maintenance which is consistent with our hypothesis for RKD-induced tillering
Cytoscape network construction
We found a number of candidates that could explain our tillering phenotype but we felt we were missing an overall picture of how indTaRKD expression was leading to our phenotype of interest. To address this we decided to construct a coexpression network.
Due to poor functional annotation in wheat and lacking an easily accessible coexpression table we decided to construct our RKD network in rice using the OsRKD DEGs as a template. Below is a brief description of the analysis and how the coexpression network was constructed.
In order to construct a coexpression network, we took a high confidence coexpression table produced by the netminer algorithm. From we only took the interactions that had an OsRKD DEG (lfcThreshold=1, p.value = 0.05) with over 0.15 weighted coexpression correlation. From this we loaded the table into cytoscape and layed it out using a prefuse force directed layout. We manually curated clusters and did enrichment tests for GO terms and 1000bp promoter motifs. An image is presented below but additional network options can be visualised using the search terms below in nDEX.
Figure 10A. Normalised expression of ERF113 (DEG) and its homeolog ERF115 across samples (non-DEG)
DEGs network options - “OsRKD OsRKD DEGs”, “OsRKD TaRKD DEGs”, “OsRKD Microspore DEGs”, “OsRKD Tillering DEGs”, “OsRKD NGR5 DEGs”
GeneNet network options - “OsRKD TaRKD geneNet 0.15”, “OsRKD Microspore geneNet 0.15”, “OsRKD Tillering geneNet 0.15”, “OsRKD NGR5 geneNet 0.15”
Tissue expression options - “OsRKD Pre-emergence inflorescence”, “OsRKD Post-emergence inflorescence”, “OsRKD Anther”, “OsRKD Pistil”, “OsRKD Embryo 25 DAP”, “OsRKD Endosperm 25 DAP”, “OsRKD Seed 5 DAP”, “OsRKD Seed 10 DAP”, “OsRKD Leaves 20 days”, “OsRKD Seedling 4 leaf stage”, “OsRKD Shoots 25 DAP”
Functional term enrichment results for clusters
Photoinhibition cluster
De-novo: MYBS2 and methionine biosynthetic genes
Known: MYB related and ERF
Jasmonate cluster
De-novo: CAMTA and RRM
Known: CAMTA, WRKY and RKD2
In general very high number of motifs compared to other cluster. This suggest is the main transcriptional regulatory cluster of the network
Signal transduction
De-novo: none
Known: E2F, CAMTA and MADS
Brassinosteroid
De-novo: Highly signficant squamosa like promoter (SPL11)
Known: MYB and P53 like
Final note
There seems to be two major processes being detected in the transcriptional networks related to RKD expression in the grasses. The first regulates transcription and cell division which is consistent with RKDs role as a regulator of early embryonic development and would be dependent on MYB and Sigma70 factors implicated in the network. The second may explain RKD induced tillering and seems to be dependent on AP2-ERF factors and how they influence phytohormone biosynthesis that increases formation of axillary meristems. This is particularly related to the phytohormone jasmonate which has been shown to play a role in regeneration and response to wounding.